• Коэффициенты сходства и различия
• Неметрическое многомерное шкалирование
• Иерархическая кластеризация

Вы сможете

• Выбирать подходящий для данных коэффициент сходства/различия
• Представлять многомерные данные в меньшем числе измерений при помощи неметрического многомерного шкалирования
• Строить дендрограммы при помощи подходящего метода аггрегации

• Используются в качестве исходных данных для многих видов многомерных анализов, в т.ч. для неметрического многомерного шкалирования и некоторых видов кластерного анализа
• Из сходств можно получить расстояния и наоборот
• Свои коэффициенты для количественных и качественных признаков

Метрики без стандартизации

• Евклидово расстояние

Неевклидовы метрики

• Квадрат Евклидова расстояния

• Манхэттеновское расстояние

Можно стандартизовать исходные данные

Простая стандартизация не подходит (счет, не может быть среднее 0)
Можно использовать трансформации: - корень, корень 4-й степени - логарифмирование со сдвигом (log10(1+n))

Можно использовать коэффициенты со стандартизацией

• Канберра (метрика) \(d = {\sum \frac {|x _{ik} - x _{jk}|} {|x _{ik}|+|x _{jk}|}}\)

• хи-квадрат (метрика) \(\chi^2 = \sqrt {\sum{ {\frac {1} {c _k}} (x _{ik} - x _{jk})^2 }}\)

• хи-квадрат (метрика) \(\chi^2 = \sqrt {\sum{ {\frac {1} {c _k}} (x _{ik} - x _{jk})^2 }}\)

• коэффициент Брея-Куртиса (не метрика) \(BC _{ij} = \frac { 2C _{ij}} {S _i + S _j}\)

• Евклидово расстояние \(d = \sqrt {\sum{(x _{ik} + x _{jk})^2}}\)

Коэффициенты для смешанных данных

• расстояние Говера

library(DAAG)
data(possum)
colnames(possum)

##  [1] "case"     "site"     "Pop"      "sex"      "age"      "hdlngth" 
##  [7] "skullw"   "totlngth" "taill"    "footlgth" "earconch" "eye"     
## [13] "chest"    "belly"

sum(is.na(possum))

## [1] 3

possum[!complete.cases(possum), ]

##      case site Pop sex age hdlngth skullw totlngth taill footlgth
## BB36   41    2 Vic   f   5    88.4   57.0       83  36.5       NA
## BB41   44    2 Vic   m  NA    85.1   51.5       76  35.5     70.3
## BB45   46    2 Vic   m  NA    91.4   54.4       84  35.0     72.8
##      earconch  eye chest belly
## BB36     40.3 15.9  27.0  30.5
## BB41     52.6 14.4  23.0  27.0
## BB45     51.2 14.4  24.5  35.0

# Добавим названия сайтов
possum$site <- factor(possum$site, levels = 1:7, labels = c("Cambarville", 
    "Bellbird", "Whian Whian", "Byrangery", "Conondale ", "Allyn River", 
    "Bulburin"))

Отберем переменные, с которыми будем работать

colnames(possum)

##  [1] "case"     "site"     "Pop"      "sex"      "age"      "hdlngth" 
##  [7] "skullw"   "totlngth" "taill"    "footlgth" "earconch" "eye"     
## [13] "chest"    "belly"

possumc <- possum[complete.cases(possum), c(3:4, 5:14)]

Построим ординацию поссумов на основе их сходства по морфометрии и возрасту.

Функция metaMDS трансформирует и стандартизует данные, а затем итеративно подбирает координаты поссумов в новом пространстве (двумерном по умолчанию).

library(vegan)
ord_euclid <- metaMDS(possumc[, 3:10], distance = "euclid")

## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.0382 
## Run 1 stress 0.0403 
## Run 2 stress 0.0403 
## Run 3 stress 0.0441 
## Run 4 stress 0.0382 
## ... procrustes: rmse 0.000539  max resid 0.00428 
## *** Solution reached

stress - оценивает, насколько были искажены исходные расстояния между объектами при снижении размерности

ord_euclid$stress

## [1] 0.0382

• Эмпирическое правило предложено Краскалом
• S > 20% плохо
• S = 10% нормально
• S < 5% хорошо
• S = 0 прекрасно

При помощи ggplot2 постройте график неметрического многомерного шкалирования.

Для графика используйте координаты точек ord_euclid$points и исходные данные.

Раскрасьте график по значениям переменных Pop и age

Изобразите поссумов разного пола на разных панелях

library(ggplot2)
theme_set(theme_bw(base_size = 12) + theme(legend.key = element_blank()))
ord_euclid_points <- data.frame(ord_euclid$points, possumc)
gg <- ggplot(ord_euclid_points, aes(x = MDS1, y = MDS2)) + geom_point() + 
    facet_wrap(~sex) + theme(legend.key = element_blank())

Примеры графиков nMDS ординации

gg + aes(colour = Pop)

gg + aes(colour = age)

Постройте nMDS ординацию при помощи евклидова расстояния, без стандартизации

Воспользуйтесь справкой к функции metaMDS(), чтобы узнать, какие аргументы потребуется изменить.

Стресс сильно вырос, эта ординация хуже. Почему так произошло?

ord_raw <- metaMDS(possumc[, 3:10], dis = "euclide", autotransform = FALSE)

## Run 0 stress 0.101 
## Run 1 stress 0.101 
## ... New best solution
## ... procrustes: rmse 0.000036  max resid 0.000253 
## *** Solution reached

ord_raw$stress

## [1] 0.101

Для сравнения, стресс был

ord_euclid$stress

## [1] 0.0382

Прокрустово преобразование

procrustes(что_стало, что_было)

proc <- procrustes(ord_raw, ord_euclid)
proc

## 
## Call:
## procrustes(X = ord_raw, Y = ord_euclid) 
## 
## Procrustes sum of squares:
## 5.67e+03

plot(proc)

ord_euclid_points <- data.frame(ord_euclid$points, possumc)
grid.arrange(gg %+% ord_euclid_points + aes(colour = Pop), gg %+% ord_euclid_points + 
    aes(colour = age), ncol = 2)

Морфометрия самок поссумов

library(DAAG)
data(fossum)
# создадим 'говорящие' имена строк
rownames(fossum) <- paste(fossum$Pop, rownames(fossum), sep = "_")
fossumc <- fossum[complete.cases(fossum), 5:14]

Методы класстеризации на основании расстояний (о них сегодня)

• Метод ближайшего соседа
• Метод отдаленного соседа
• Метод среднегруппового расстояния
• Метод Варда
• и т.д. и т.п.

Методы кластеризации на основании признаков

• Метод максимальной бережливости
• Метод максимального правдоподобия

Результат кластеризации зависит от

• коэффициента сходства-различия
• от алгоритма кластеризации

В начале лекции мы обнаружили, что евклидово расстояние, расчитанное по стандартизованным данным, хорошо разделяет поссумов.

d <- dist(x = scale(fossumc), method = "euclidean")

Давайте построим деревья при помощи нескольких алгоритмов кластеризации и сравним их.

Особенности:

• Может быть сложно интерпретировать, если нужны группы
• объекты на дендрограмме часто не образуют четко разделенных групп
• часто получаются цепочки кластеров (объекты присоединяются как бы по-одному)
• Хорош для выявления градиентов

hc_single <- hclust(d, method = "single")
library(ape)
ph_single <- as.phylo(hc_single)
plot(ph_single, type = "phylogram", cex = 0.7)
axisPhylo()

Особенности:

• На дендрограмме образуется много отдельных некрупных групп
• Хорош для поиска дискретных групп в данных

ph_compl <- as.phylo(hclust(d, method = "complete"))
plot(ph_compl, type = "phylogram", cex = 0.7)
axisPhylo()

Особенности:

• UPGMA и WUPGMС иногда могут приводить к инверсиям на дендрограммах

ph_avg <- as.phylo(hclust(d, method = "average"))
plot(ph_avg, type = "phylogram", cex = 0.7)
axisPhylo()

Особенности:

• метод годится и для неевклидовых расстояний несмотря на то, что внутригрупповая дисперсия расстояний рассчитывается так, как будто это евклидовы расстояния

ph_w2 <- as.phylo(hclust(d, method = "ward.D2"))
plot(ph_w2, type = "phylogram", cex = 0.7)
axisPhylo()

Кофенетическое расстояние - расстояние между объектами на дендрограмме

Кофенетическую корреляцию можно рассчитать как пирсоновскую корреляцию (обычную) между матрицами исходных и кофенетических расстояний между всеми парами объектов

Метод, который дает наибольшую кофенетическую корреляцию дает кластеры лучше всего отражающие исходные данные

c_single <- cophenetic(ph_single)
c_compl <- cophenetic(ph_compl)
c_avg <- cophenetic(ph_avg)
c_w2 <- cophenetic(ph_w2)
cor(d, as.dist(c_single))

## [1] 0.679

cor(d, as.dist(c_compl))

## [1] 0.524

cor(d, as.dist(c_avg))  # лучше всех отражает структуру данных

## [1] 0.742

cor(d, as.dist(c_w2))

## [1] 0.541

• Можно субъективно, на любом выбранном уровне. Главное, чтобы кластеры были осмысленными и интерпретируемыми.
• Можно выбрать, глядя на распределение расстояний ветвления
• Можно оценить вероятность разделения на кластеры при помощи бутстрепа

library(pvclust)

# итераций должно быть 1000 и больше здесь мало для скорости
set.seed(42)
cl_boot <- pvclust(scale(t(fossumc)), method.hclust = "average", nboot = 50, 
    method.dist = "euclidean")

## Bootstrap (r = 0.5)... Done.
## Bootstrap (r = 0.6)... Done.
## Bootstrap (r = 0.7)... Done.
## Bootstrap (r = 0.8)... Done.
## Bootstrap (r = 0.9)... Done.
## Bootstrap (r = 1.0)... Done.
## Bootstrap (r = 1.1)... Done.
## Bootstrap (r = 1.2)... Done.
## Bootstrap (r = 1.3)... Done.
## Bootstrap (r = 1.4)... Done.

plot(cl_boot)
pvrect(cl_boot)

webpage <- "http://evolution.genetics.washington.edu/book/primates.dna"
primates.dna <- read.dna(webpage)
d_pri <- dist.dna(primates.dna)
hc_pri <- hclust(d_pri, method = "average")
ph_pri <- as.phylo(hc_pri)
plot(ph_pri)
axisPhylo()

• Неметрическое многомерное шкалирование - способ снижения размерности, сохраняющий ранги расстояний между объектами
• Направления на графике многомерного шкалирования можно интерпретировать произвольным образом в зависимости от изменения других переменных (не обязательно вдоль осей)
• Результат многомерного шкалирования зависит от выбора коэффициента различия

• Borcard, D., Gillet, F., Legendre, P., 2011. Numerical ecology with R. Springer.
• Legendre, P., Legendre, L., 2012. Numerical ecology. Elsevier.
• Oksanen, J., 2011. Multivariate analysis of ecological communities in R: vegan tutorial. R package version 2–0.
• Quinn, G.G.P., Keough, M.J., 2002. Experimental design and data analysis for biologists. Cambridge University Press.

Borcard, D., Gillet, F., Legendre, P., 2011. Numerical ecology with R. Springer.

Как работает UPGMA можно посмотреть здесь: - http://www.southampton.ac.uk/~re1u06/teaching/upgma/

• pvclust: An R package for hierarchical clustering with p-values [WWW Document], n.d. URL http://www.sigmath.es.osaka-u.ac.jp/shimo-lab/prog/pvclust/ (accessed 11.7.14).

Для анализа молекулярных данных: - Paradis, E., 2011. Analysis of Phylogenetics and Evolution with R. Springer.